%
% vcycle_block.m
%
%   recursive vcycle code using block smoother
%
function v = vcycle_block(u,f,level);
  
  %disp(sprintf('level = %5i',level));
  
  % data structure with the matricies, and parameters
  %
  global mg_matrix
  global mg_param

  % multigrid parameters
  %
  nu1      = mg_param.nu1;
  nu2      = mg_param.nu2;
  minlevel = mg_param.minlevel;
  
  % block smoother parameters
  %
  bx = mg_param.bx;
  by = mg_param.by;
  sx = mg_param.sx;
  sy = mg_param.sy;
  
  
  % if on the last level solve the muthafucka 
  %
  if( level == minlevel )
      B = mg_matrix.A{level};
      if (condest(B)>1.e12) % this is for pressure or phi
          sz = size(B);
          B = [B; ones(1,sz(2))];
          f = [f;0];
      end
      v = B\f;      
      %v = mg_matrix.A{level}\f;    
  else
    
    % presmooth
    %
    Nx = 2^level;  % size of the current grid
    Ny = Nx;
    v = blocksmooth(u,mg_matrix.A{level},f,nu1,Nx,Ny,bx,by,sx,sy);
        
    % compute residual
    %
    r = f - mg_matrix.A{level}*v;
    
    % restrict
    %
    fc = mg_matrix.Ir{level}*r;
    
    % initialize
    %
    uc = 0.0*fc;
    
    % solve the error equation, recursively
    %
    ec = vcycle_block(uc,fc,level-1);
    
    % correct
    %
    v = v + mg_matrix.Ip{level} * ec;
  
    % post smoothing
    %
    v = blocksmooth(v,mg_matrix.A{level},f,nu2,Nx,Ny,bx,by,sx,sy);  
  end
  
  
